#include "field.h"

double care(double a)
{
	return(a*a); 
}

void field_nonlinear_1(Field       *f_rx,
					   const Field *f_source1, 
					   const Field *f_source2,
					   const Field *f_rho,
					   const double inv_dx, 
					   const double inv_dy
					 )
{
	CHECK_FIELD(f_rx);
	CHECK_FIELD(f_source1); 
	CHECK_FIELD(f_source2); 
	CHECK_FIELD(f_rho);
	CHECK_FIELDS(f_source1, f_source2);
	CHECK_FIELDS(f_source1, f_rho); 
	CHECK_FIELDS(f_rx, f_rho); 
	{
		const  int        nx = f_rx->nx; 
		const  int        ny = f_rx->ny; 
		const  int       nx1 = nx-1; 
		const  int       ny1 = ny-1; 
		double **        rx  = ARRAY_OF(f_rx); 
		const double   **rho = CONST_ARRAY_OF(f_rho);
		const double   **rhu = CONST_ARRAY_OF(f_source1); 
		const double   **rhv = CONST_ARRAY_OF(f_source2); 
		
		int i,j; 
		
		for(j=0;j<ny;++j)
        {
            const int jm = (j <= 0)   ? ny1 : j-1;
            const int jp = (j >= ny1) ? 0   : j+1;
            
            for(i=0;i<nx;++i)
            {
                const int im = (i <= 0)   ? nx1 : i-1;
                const int ip = (i >= nx1) ? 0   : i+1;
				
				rx[j][i] = inv_dx*(care(0.5*(rhu[j][ip]+rhu[j][i]))/rho[j][ip]-care(0.5*(rhu[j][i]+rhu[j][im]))/rho[j][i])+inv_dy*((0.25*(rhu[jp][i]+rhu[j][i])*(rhv[j][i]+rhv[j][ip]))/(0.25*(rho[jp][ip]+rho[jp][i]+rho[j][i]+rho[j][ip]))-(0.25*(rhu[j][i]+rhu[jm][i])*(rhv[jm][i]+rhv[jm][ip]))/(0.25*(rho[j][i]+rho[j][ip]+rho[jm][ip]+rho[jm][i]))); 
			}
		}
		
	}
	
}
	
void field_nonlinear_2(Field       *f_ry,
					   const Field *f_source1, 
					   const Field *f_source2,
					   const Field *f_rho,
					   const double inv_dx, 
					   const double inv_dy
					   )
{
	CHECK_FIELD(f_ry); 
	CHECK_FIELD(f_source1); 
	CHECK_FIELD(f_source2); 
	CHECK_FIELD(f_rho);
	CHECK_FIELDS(f_source1, f_source2);
	CHECK_FIELDS(f_source1, f_rho); 
	CHECK_FIELDS(f_ry, f_rho); 
	{
		const int         nx = f_ry->nx; 
		const int         ny = f_ry->ny; 
		const int        nx1 = nx-1; 
		const int        ny1 = ny-1; 
		double **         ry = ARRAY_OF(f_ry); 
		const double   **rho = CONST_ARRAY_OF(f_rho);
		const double   **rhu = CONST_ARRAY_OF(f_source1); 
		const double   **rhv = CONST_ARRAY_OF(f_source2); 
		
		int i,j; 
		
		for(j=0;j<ny;++j)
        {
            const int jm = (j <= 0)   ? ny1 : j-1;
            const int jp = (j >= ny1) ? 0   : j+1;
            
            for(i=0;i<nx;++i)
            {
                const int im = (i <= 0)   ? nx1 : i-1;
                const int ip = (i >= nx1) ? 0   : i+1;
				
				ry[j][i] = inv_dx*((0.25*(rhv[j][i]+rhv[j][ip])*(rhu[j][i]+rhu[jp][i]))/(0.25*(rho[jp][ip]+rho[jp][i]+rho[j][i]+rho[j][ip]))-((0.25*(rhv[j][i]+rhv[j][im])*(rhu[jp][im]+rhu[j][im]))/(0.25*(rho[jp][im]+rho[j][im]+rho[j][i]+rho[jp][i]))))+inv_dy*(care(0.5*(rhv[jp][i]+rhv[j][i]))/rho[jp][i]-care(0.5*(rhv[jm][i]+rhv[j][i]))/rho[j][i]); 
			}
		}
		
	}
	
}
     